% Function:
function [I, R, RR] = Backward(Y, X, l)
    if nargin ~= 3
        error('Usage: Backward(Y, X, l)');
    end

    [m, n] = size(X);   % m: number of samples, n: number of features

    % initialize collection
    F  = [ones(m, 1), X];   % the training set
    I  = zeros(1, n - l);   % index of discarded features
    R  = zeros(1, n - l);   % average residuals
    RR = zeros(1, n - l);   % R^2 statistic

    idx = 1 : n;

    k = n;
    while k > l
        r   = zeros(n);
        rr  = ones(n);

        % search feature i to n
        for i = 1 : n
            % feature i is already discarded
            if any(i == I)
                continue;
            end

            slct = idx;
            slct(slct == i) = [];
            slct = slct + 1;
            slct = [1, slct];

            [~, ~, t, ~, stats] = regress(Y, F(:, slct), 0.05);
            r(i)  = mean(t .^ 2);
            rr(i) = stats(1);
        end

        % select the feature that has the max R^2 statistic
        mini = find(rr == min(rr));
        mini = mini(1);

        % udpate collection
        i     = n - k + 1;
        I(i)  = mini;
        R(i)  = r(mini);
        RR(i) = rr(mini);

        idx(idx == mini) = [];

        % update loop variable
        k = k - 1;
    end
end
